Deconvolution by FIR.rst

%pylab inline
import numpy as np
import scipy.signal as dsp
from palettable.colorbrewer.qualitative import Dark2_8

colors = Dark2_8.mpl_colors
rst = np.random.RandomState(1)

Populating the interactive namespace from numpy and matplotlib

Design of a digital deconvolution filter (FIR type)

from PyDynamic.model_estimation.fit_filter import LSFIR
from PyDynamic.misc.SecondOrderSystem import *
from PyDynamic.misc.testsignals import shocklikeGaussian
from PyDynamic.misc.filterstuff import kaiser_lowpass, db
from PyDynamic.uncertainty.propagate_filter import FIRuncFilter
from PyDynamic.misc.tools import make_semiposdef

# parameters of simulated measurement
Fs = 500e3
Ts = 1 / Fs

# sensor/measurement system
f0 = 36e3; uf0 = 0.01*f0
S0 = 0.4; uS0= 0.001*S0
delta = 0.01; udelta = 0.1*delta

# transform continuous system to digital filter
bc, ac = sos_phys2filter(S0,delta,f0)
b, a = dsp.bilinear(bc, ac, Fs)

# Monte Carlo for calculation of unc. assoc. with [real(H),imag(H)]
f = np.linspace(0, 120e3, 200)
Hfc = sos_FreqResp(S0, delta, f0, f)
Hf = dsp.freqz(b,a,2*np.pi*f/Fs)[1]

runs = 10000
MCS0 = S0 + rst.randn(runs)*uS0
MCd  = delta+ rst.randn(runs)*udelta
MCf0 = f0 + rst.randn(runs)*uf0
HMC = np.zeros((runs, len(f)),dtype=complex)
for k in range(runs):
    bc_,ac_ = sos_phys2filter(MCS0[k], MCd[k], MCf0[k])
    b_,a_ = dsp.bilinear(bc_,ac_,Fs)
    HMC[k,:] = dsp.freqz(b_,a_,2*np.pi*f/Fs)[1]

H = np.r_[np.real(Hf), np.imag(Hf)]
uAbs = np.std(np.abs(HMC),axis=0)
uPhas= np.std(np.angle(HMC),axis=0)
UH= np.cov(np.hstack((np.real(HMC),np.imag(HMC))),rowvar=0)
UH= make_semiposdef(UH)

Problem description

Assume information about a linear time-invariant (LTI) measurement system to be available in terms of its frequency response values H(jω) at a set of frequencies together with associated uncertainties:

u(H) = (u(|H(jω₁)|), …, u(|H(jω_N)|), u(∠H(jω₁)), …, u(∠H(jω_N)))

figure(figsize=(16,8))
errorbar(f*1e-3, np.abs(Hf), uAbs, fmt=".", color=colors[0])
title("measured amplitude spectrum with associated uncertainties")
xlim(0,50)
xlabel("frequency / kHz",fontsize=20)
ylabel("magnitude / au",fontsize=20);

figure(figsize=(16,8))
errorbar(f*1e-3, np.angle(Hf), uPhas, fmt=".", color=colors[1])
title("measured phase spectrum with associated uncertainties")
xlim(0,50)
xlabel("frequency / kHz",fontsize=20)
ylabel("phase / rad",fontsize=20);

Simulated measurement

Measurements with this system are then modeled as a convolution of the system's impulse response

h(t) = ℱ^− 1(H(jω))

with the input signal x(t), after an analogue-to-digital conversion producing the measured signal

y[n] = (h * x)(t_n)  n = 1, …, M

# simulate input and output signals
time = np.arange(0, 4e-3 - Ts, Ts)
#x = shocklikeGaussian(time, t0 = 2e-3, sigma = 1e-5, m0=0.8)
m0 = 0.8; sigma = 1e-5; t0 = 2e-3
x = -m0*(time-t0)/sigma * np.exp(0.5)*np.exp(-(time-t0) ** 2 / (2 * sigma ** 2))
y = dsp.lfilter(b, a, x)
noise = 1e-3
yn = y + rst.randn(np.size(y)) * noise

figure(figsize=(16,8))
plot(time*1e3, x, label="system input signal", color=colors[0])
plot(time*1e3, yn,label="measured output signal", color=colors[1])
legend(fontsize=20)
xlim(1.8,4); ylim(-1,1)
xlabel("time / ms",fontsize=20)
ylabel(r"signal amplitude / $m/s^2$",fontsize=20);

Design of the deconvolution filter

The aim is to derive a digital filter with finite impulse response (FIR)

$$g(z) = \sum_{k=0}^K b_k z^{-k}$$

such that the filtered signal

x̂[n] = (g * y)[n]  n = 1, …, M

is an estimate of the system's input signal at the discrete time points.

Publication

• Elster and Link "Uncertainty evaluation for dynamic measurements modelled by a linear time-invariant system" Metrologia, 2008
• Vuerinckx R, Rolain Y, Schoukens J and Pintelon R "Design of stable IIR filters in the complex domain by automatic delay selection" IEEE Trans. Signal Process. 44 2339–44, 1996

Determine FIR filter coefficients such that

H(jω)g(e^jω/F_s) ≈ e^{− jωn₀/F_s}  for  |ω| ≤ ω₁

with a pre-defined time delay n₀ to improve the fit quality (typically half the filter order).

Consider as least-squares problem

(y − Xb)^TW^− 1(y − Xb)

with - y real and imaginary parts of the reciprocal and phase shifted measured frequency response values - X the model matrix with entries e^{− jkω/Fs} - b the sought FIR filter coefficients - W a weighting matrix (usually derived from the uncertainties associated with the frequency response measurements

Filter coefficients and associated uncertainties are thus obtained as

b = (X^TW^− 1X)^− 1X^TW^− 1y

u_b = (X^TW^− 1X)^− 1X^TW^− 1U_yW^− 1X(X^TW^− 1X)^− 1

# Calculation of FIR deconvolution filter and its assoc. unc.
N = 12; tau = N//2
bF, UbF = LSFIR(H,N,tau,f,Fs,UH=UH)

Least-squares fit of an order 12 digital FIR filter to the reciprocal of a frequency response given by 400 values and propagation of associated uncertainties. Final rms error = 1.545423e+01

figure(figsize=(16,8))
errorbar(range(N+1), bF, np.sqrt(np.diag(UbF)), fmt="o", color=colors[3])
xlabel("FIR coefficient index", fontsize=20)
ylabel("FIR coefficient value", fontsize=20);

In order to render the ill-posed estimation problem stable, the FIR inverse filter is accompanied with an FIR low-pass filter.

Application of the deconvolution filter for input estimation is then carried out as

x̂[n − n₀] = (g * (g_low * y)[n]

with point-wise associated uncertainties calculated as

u²(x̂[n − n₀] = b^TU_xlow[n]b + x_low^T[n]U_bx_low[n] + trace(U_xlow[n]U_b)

fcut = f0+10e3; low_order = 100
blow, lshift = kaiser_lowpass(low_order, fcut, Fs)
shift = -tau - lshift

figure(figsize=(16,10))
HbF = dsp.freqz(bF,1,2*np.pi*f/Fs)[1]*dsp.freqz(blow,1,2*np.pi*f/Fs)[1]
semilogy(f*1e-3, np.abs(Hf), label="measured frequency response")
semilogy(f*1e-3, np.abs(HbF),label="inverse filter")
semilogy(f*1e-3, np.abs(Hf*HbF), label="compensation result")
legend();

xhat,Uxhat = FIRuncFilter(yn,noise,bF,UbF,shift,blow)

figure(figsize=(16,8))
plot(time*1e3,x, label='input signal')
plot(time*1e3,yn,label='output signal')
plot(time*1e3,xhat,label='estimate of input')
legend(fontsize=20)
xlabel('time / ms',fontsize=22)
ylabel('signal amplitude / au',fontsize=22)
tick_params(which="both",labelsize=16)
xlim(1.9,2.4); ylim(-1,1);

figure(figsize=(16,10))
plot(time*1e3,Uxhat)
xlabel('time / ms',fontsize=22)
ylabel('signal uncertainty / au',fontsize=22)
subplots_adjust(left=0.15,right=0.95)
tick_params(which='both', labelsize=16)
xlim(1.9,2.4);

Basic workflow in PyDynamic

Fit an FIR filter to the reciprocal of the measured frequency response

from PyDynamic.model_estimation.fit_filter import LSFIR
bF, UbF = LSFIR(H,N,tau,f,Fs,verbose=False,UH=UH)

with

• H the measured frequency response values
• UH the covariance (i.e. uncertainty) associated with real and imaginary parts of H
• N the filter order
• tau the filter delay in samples
• f the vector of frequencies at which H is given
• Fs the sampling frequency for the digital FIR filter

Propagate the uncertainty associated with the measurement noise and the FIR filter through the deconvolution process

xhat,Uxhat = FIRuncFilter(yn,noise,bF,UbF,shift,blow)

with

• yn the noisy measurement
• noise the std of the noise
• shift the total delay of the FIR filter and the low-pass filter
• blow the coefficients of the FIR low-pass filter